import math
import matplotlib.pyplot as plt
import numpy as np
import sympy as sy
from sympy import *
def get_temp(a,b,n):
    temp = (b-a) /n
    return temp
#-------------相关参数-------------
a = 1
b = 2
n = 8
eps = 10** -5 #精度
T = []
S = []
C = []
R = []


#构建原函数
temp = get_temp(a,b,n)
scalex = np.arange(a,b+temp,temp)
ylist = []
for i in scalex:
    if(i == 0):
        ylist.append(1);
    else:
       ylist.append(math.sin(i) / i)
def func(x):
    if(x==0):
      return 1
    else:
      return math.sin(x) / x


#复合梯形公式
def ti(x,fx,n):
   h = get_temp(a,b,n)
   res = 0
   for i in range(0,len(x)):
       if(i==0 or i == len(x)-1):
           res+=fx[i]
       else:
           res+= 2*fx[i]
   res = (h/2) * res
   return res


#复合辛普森
def xinpusen(a,b,n):
    h = (b-a) / n
    res = func(a) - func(b)
    x = a
    for i in range(1,n+1):
        x = x + (h/2)
        res += 4*(func(x))
        x = x + (h/2)
        res += 2*(func(x))
    res = (h/6)*res
    return res

    
#龙贝格
def romberg(a,b,eps,func):
    h = b-a
    #初始化T
    T.append(h * (func(a) + func(b)) /2 )
    ep = eps + 1
    m = 0       #加速M次其实是二分M次，若二分次数到达相应次数，则加速
    while(ep>=eps):
        m = m+1
        t = 0
        #每次重新计算t的值
        #每次二分都加上 2**(m-1)个点
        for i in range(2**(m-1)-1):
            #可以得到每次二分之后的点都是2的二分次数作为分母。相应的奇数作为分子，
            t = t +func(a+(2*(i+1)-1)*h / 2**m)*h/2**m
            #计算T2n
        t = t + T[-1]/2
        T.append(t)
            #计算每一次加速之后的函数
        if m>=1:
                S.append((4**m*T[-1]-T[-2])/(4**m-1))
        if m>=2:
                C.append((4**m*S[-1]-S[-2])/(4**m-1))
        if m>=3:   
                R.append((4**m*C[-1]-C[-2])/(4**m-1))
        if m>4:
                ep=abs(10*(R[-1]-R[-2]))
#自适应积分
def auto(a,b,relresult):
    mid = (b-a) /2 
    left = xinpusen(a,mid,2)
    right = xinpusen(mid,b,2)
    if(abs(left + right - relresult)<eps):
        return left + right
    else:
        return auto(a,mid,left) + auto(mid,b,right)
#复合梯形公式余项
def get_TR(a,b,n):
    temp = get_temp(a,b,n)
    RT = (1/12) * (temp**2) * (1/3)
    return RT

#复合辛普森公式余项
def get_SR(a,b,n):
    temp = get_temp(a,b,n/2)
    RS = (1/180) * (math.pow(temp/2,4)) * (1/5)
    return RS

#龙贝格公式余项
def get_LR(a,b,n):
    temp =get_temp(a,b,n/6)
    LR = (1/40320) * (1/9*(math.pow(b,9)- math.pow(a,9))  - ((b-a)/840 * (41*math.pow(a,8) + 216 * math.pow((b-a)/6,8) + 27 * math.pow((b-a)/3,8) + 
    272* math.pow((b-a)/2,8) + 27 * math.pow(2*(b-a)/3,9) + 216 * math.pow(5*(b-a)/6,8) + 41 * math.pow(b,8))))
    return LR


#复合梯形公式
res1 = ti(scalex,ylist,n)
print("复合梯形的结果为:",res1)
print("复合梯形公式的余项为:",get_TR(a,b,n))

#复合辛普森公式
res2 = xinpusen(a,b,n)
print("复合辛普森的结果为:",res2)
print("复合辛普森公式的余项为:",get_SR(a,b,n))

#龙贝格
romberg(a,b,eps,func)
print("龙贝格的结果为:",R[-1])
print("龙贝格公式的余项为:",get_LR(a,b,n))


#自适应积分
res3 = auto(a,b,0.5619649373692952)
print("自适应积分的结果:",res3)
#plt.plot(scalex,ylist,color = 'red',label = 'six/x')
#plt.legend()
#plt.show()